Identification of arches in 2D granular packings 
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SO 1 

■ We identify arches in a bed of granular disks generated by a molecular dynamic-type simulation. 

We use the history of the deposition of the particles to identify the supporting contacts of each 
particle. Then, arches are defined as sets of mutually stable disks. Different packings generated 
through tapping are analyzed. The possibility of identifying arches from the static structure of a 
r-{ ' deposited bed, without any information on the history of the deposition, is discussed. 

hj ; 

' I- INTRODUCTION 

Arches are multi-particle structures encountered in non-sequentially deposited granular beds. In some sense, arches 
are to granular packings what clusters are to colloids: they represent collective structures whose existence and behavior 
determine many properties of the system as a whole. During the settlement of a granular bed, some particles come 
. to rest simultaneously by contacting and supporting each other. These mutually stabilized sets of particles are called 
arches or bridges 1]. Of course, these arches are themselves supported by some of the surrounding particles in the 
system: the arch bases. Typically, 70% of the particles are part of arches in a granular packing of hard spheres 0]; 
and around 60% in packings of hard disks 3]. Arch formation is crucial in the jammingof granular flows driven by 
gravity (3, 0, IE an d it has also been proposed as a mechanism for size segregation H,H ■ 

Arching is the collective process associated to the appearance of voids that lower the packing fraction of the sample. 
Moreover, arching is directly related to the reduction of particle-particle contacts in the assembly, which determines 
the coordination number pi Il0|. In 2D, for example, the mean support number (z) S upport — i-e., the coordination 
' number that accounts 0Ii hL f° r the contacts that serve as support of at least one grain — can be obtained from the arch 
£>. , size distribution n(s) as 3J 



T3 



-a 
c 
o 



- 1—1 

X 



l N 

-Yn(s) 



s=l 



(1) 



o : 

oo . 

■ {-^support 

O ' 

vo : . 

Here N is the total number of particles in the packing and n(s) is the number of arches consisting of s grains, 
"j***; ■ The number of particles that do not form part of any arch corresponds to n(s = 1). This mean support number 
coincides with the mean coordination number (z) when the packing does not have any non-supporting contacts. In 
non-sequentially deposited beds — especially for soft particles — there will be contacts that are not essential to the 
I ' stability of any of the two touching particles. 

A systematic study of arching is particularly complex despite the simplicity associated to the concept of arch. 
The typical structural properties of granular materials are mainly connected with the topological complexity of the 
contact network. Concepts like "force chain" and "force propagation" have been profusely studied [Tlj , but many open 
L" ' questions still remain. The statistical properties of the arching process could give some insight into these complex 
problems. 

Arches form dynamically during non-sequential particle deposition. Given a settled granular sample, identifying 
if two particles belong to the same arch is often impossible without knowing the history of the deposition process. 
Assuming that we are able to identify (experimentally) the contacts of each particle, it is impossible to know which of 
these contacts are responsible for supporting each particle in place against gravity. For convex hard particles in 2D, 
only two of the contacts of a given particle provide its stability (in 3D this number is three) . The supporting contacts 
are the first two (three in 3D) contacts made by the given particle that provide stability against the external force 
that drives deposition (e.g. gravity). Any contact made afterwards will not provide the essential stability assuming 
that the already formed supporting contacts persist. Of course, if one removes a supporting contact, it is possible 
that a non-supporting contact of the particle becomes a supporting contact, hence some non-supporting contacts may 
provide secondary (alternative) stability to the particle. 

Finding the supporting contacts of a particle is simple in a pseudo dynamics simulation approach since these 
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contacts are the essence of the simulation algorithm pi 13. IToL Il2|. However, in a realistic granular simulation 



well as in experiments — one needs to track the particle collisions and analyze stability at every instant to catch the 
supporting contacts of every particle in the moment they first appear. 

In this work, we report a numerical study on the supporting contacts that are created during the non-sequential 
deposition of a granular bed. Then, arches are straightforwardly obtained and analyzed. We use molecular dynamics 
simulations and restrict ourselves to a 2D system of inelastic particles in order to reduce the computational cost. We 
show that our identification of supporting contacts is useful in deciding when the packing has settled, which is a major 
issue in these type of simulations. We compare our results with those obtained by pseudo dynamics methods. We also 
show that using simple criteria it is possible to identify arches to some degree without knowing the deposition history. 
This is of particular interest for experimental studies where tracking particles at high velocities is rather expensive. 

II. THE 2D SOFT-PARTICLE MOLECULAR DYNAMICS APPROACH 

We simulate the process of filling a box with grains by using a soft-particle two-dimensional molecular dynamics 
(MD). In this case, the particles (monosized disk) are subjected to the action of gravity. Particle-particle interactions 
are controlled by the particle-particle overlap £ = d — \rij\ and the velocities r^, LUi and ujj. Here, represents the 
center-to-center vector between particles i and j, d is the particle diameter and co is the particle angular velocity. 
These forces are introduced in the Newton's translational and rotational equations of motion and then numerically 
integrated by standard methods |13| . 

The contact interactions involve a normal force F n and a tangential force Ft. In order to simplify the calculus we 
use a normal force which involves a linear (Hookean) interaction between particles. 



where 



F n = k n £ - j n vf ;j 
F t = -mm(/i|F n |,| J F s |)-sign(0 

F s = -k s C - 7s < 3 
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C(t)= fvl 3 {t')dt' (5) 
Jt 

v tj =*ij ■s+^d(uj i +uj j ) (6) 

The first term in Eq. @ corresponds to a restoring force proportional to the superposition £ of the interacting 
disks and the stiffness constant k n . The second term accounts for the dissipation of energy during the contact and 
is proportional to the normal component vfj of the relative velocity pg of the disks. The restitution coefficient is 
an exponentially decaying function of the dissipation coefficient 7„ ,13], i.e., an increase in the dissipation coefficient 
leads to a nonlinear decline of the restitution coefficient. 

Equation @ provides the magnitude of the force in the tangential direction. It implements the Coulomb's criterion 
with an effective friction following a rule that selects between static or dynamic friction. Dynamic friction is accounted 
for by the friction coefficient /i. The static friction force F s [see Eq. Q] has an elastic term proportional to the relative 
shear displacement £ and a dissipative term proportional to the tangential component v\ j of the relative velocity. In 
Eq. 0, s is a unit vector normal to r.y. The elastic and dissipative contributions are characterized by k s and 7 S 
respectively. The shear displacement £ is calculated through Eq. by integrating v\ ^ from the beginning of the 
contact (i.e., t = to). The tangential interaction behaves like a damped spring which is formed whenever two grains 
come into contact and is removed when the contact finishes . 

The model presented above includes all features that turned out to be essential in the work here reported. If only 
dynamic friction forces are used in the tangential direction, these keep changing slightly but continuously once the 
disks have been deposited, making it impossible to decide whether the deposition process has definitively finished. 
The addition of static friction forces allows the deposit to reach a stable configuration. Furthermore, apart from the 
energy dissipation on normal contact, the tangential component is also responsible for energy dissipation through 
the second term in Eq. . This leads to a fast decay of the total kinetic energy of the system which final value is 
negligible in comparison with simulations without tangential dissipation. 
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III. ARCH IDENTIFICATION 

To identify arches one needs first to identify the two supporting particles of each disk in the packing. Then, arches 
can be identified in the usual way 0, 0, llOj : we first find all mutually stable particles — which we define as directly 
connected — and then we find the arches as chains of connected particles. Two disks i and j are mutually stable if i 
supports j and j supports i. 

Regardless of the way grains are deposited, they meet with and split up from other grains until they come to rest. 
While this process takes place, each time a new contact happens we can check if any of the two touching particles 
may serve to the stability of the other and save this information. Equally, we can update this data when touching 
grains come apart. Eventually, when all particles have settled, we should find that every particle has two supporting 
contacts and that no updates take place any longer. The algorithm that we follow to detect and update supporting 
contacts is described below. 

We consider that two particles in contact with particle i may provide support to it only if the two particles are one 
at each side of i and the center of mass (cm) of particle i is above the segment that joins the two contacts. We call 
contact chord to this segment. For a particle and a wall to support particle i we just need to take into account that 
the contact between the wall and i has the same y-coordinate as i. A particle in contact with the base is stable per 
se. 

If at any time a particle i has a single contact, we consider this contact as a potential first stabilizing contact 
(awaiting for the second stabilizing contact to occur) only if i has a higher y-coordinate than the contacting particle. 
A contacting wall cannot be considered a first stabilizing contact. 

We construct two arrays [nR(i) and nL(i) with i — 0, ...N — I] that store the indices of the right and left supporting 
particle of all particles in the system. If a particle has an undefined support the corresponding position in the array 
is set to —2. If one of the supporting contacts is a container's wall (or base) the corresponding position in the array 
is set to — 1. Initially, all elements of nR and nL are set to —2. After an update of all particle positions in the 
simulation, we check the status of nR and nL, and update them according to the protocol described below. The 
protocol is repeated until no changes are induced in nR and nL before a new simulation time step is advanced. 

There are six different situations that may arise depending on the stability status of a particle previous to the 
position update. Within each of these cases there are different actions to take depending on the new positions of the 
particles after the position update. A summary of the most relevant situations is shown in Fig. ^ Not all plausible 
situations are considered here because some do not occur in practice or are extremely rare. 

Case (a): Here, particle i was resting on the container's base before the position update. We simply check if the 
particle is still in contact with the base; if not, nR(i) and nL(i) are set to —2. 

Case (b): In this case, particle i was supported (before the position update took place) by other two particles which 
indices are stored in nR(i) and nL(i). We have to check if these particles are still supporting i; if not, the supporting 
particles must be redefined according to these possible situations: 

(1) nR{i) and nL(i) are still in contact with i, nR(i) is still to the right of i and nL(i) is still to the left of i, and 
the cm of i is above the contact chord — > leave nR(i) and nL[i) unchanged. 

(2) neither nR(i) nor nL(i) are in contact with i — > set nR(i) and nL(i) to —2. 

(3) either nR(i) or nL(i) is no longer in contact with i — > set the lost contact to —2 and check that the remaining 
contacting particle is still on its side (right or left) and below particle i: if not, set also this contact to —2. 

(4) nR(i) and nL(i) are still in contact with i, but either nR(i) is not to the right of i or nL(i) is not to the left 
of i — » say nR(i) is not to the right of i; then (i) set nL(i) to the contacting particle \nR(i) or nL(i)] with lower 
y-coordinate, (ii) set nR(i) — —2. 

(5) nR(i) and nL(i) are still in contact with i, nR(i) is still to the right of i and nL(i) is still to the left of i, but 
the cm of i is not above the contact chord — > set nR(i) or nL(i) to —2, the one with highest y-coordinate. 

Case (c): This case corresponds to a particle with a first potentially stabilizing contact. Let us assume that this 
corresponds to nR(i). We have to check that nR(i) is still in place and if any new contact can complete the stability 
condition. The following situations may arise: 

(1) nR(i) is still in contact with i, nR(i) is still to the right of i — > look for any particle j (or wall) contacting i on 
its left side such that the cm of i is above the contact chord; if any is found, then set nL(i) — j. 

(2) nR(i) is still in contact with i, but nR(i) is not to the right of i — > set nL(i) — nR(i) and nR(i) = —2. 

(3) nR(i) is not in contact with i — > set nR(i) = —2. 

Case (d): This case is similar to case (b) but one of the two container's walls takes the place of one of the 
supporting particles. Let us assume that nR(i) = — 1, i.e. we are considering the right wall. Again, we have to 
redefined supporting contacts according to these possible situations: 

(1) nR(i) and nL(i) are still in contact with i, nL(i) is still to the left of i, and the cm of i is above the contact 
chord — > leave nR(i) and nL(i) unchanged. 

(2) neither nR(i) nor nL(i) are in contact with i — > set nR(i) and nL(i) to —2. 
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case (a) 

nR(i) = -1 
nL(i) = -1 
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case (b) 

nR(i) = k 
nL(i) =j 
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case (c) 

nR(i) = k 
nL(i) = -2 
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nR{i) = -\ 
nL(i) = -\ 



nR(i) = -2 
nL(i) = -2 



nR(i) = k 
nL(i) =j 



nR(i) = -2 
nL(i) = -2 



nR(i) = -2 
nL{i) =j 



nR(i) = -2 
nL(i) = -2 



nR(f) = -2 
nL(i) = -2 



nR(i) = -2 
nL(i) = k 



nR(i) = -2 
nL(i) =j 



nR(i) = k; 

if i contacts 
j or wall and 
the cm is 
below the 
contact chord 
then set 
nL(i) = j or -1 
accordingly 



nR(i) = -2 
nL(i) = k 



nR(i) = -2 
nL(i) = -2 



case (d) 

nR(i) = -\ 
nL(i) =j 
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case (e) 

nRii) = -2 
nlii) = -2 
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nRii) = -2 
nL(i) =j 



nR(f) = -2 
nLii) = -2 



nR(i) = -2 
nL{f) = -2 



nR(i) = j 
nL(i) = -2 



nR(i) = -2 
riL(i) = -2 



nR(i) = -\ 
nUi) = -\ 



nR{i) = -2 
hZ(z') =j 



nR(i) = -2 
nL(i) = -2 



FIG. 1: Schematic representation of the possible states of stability that a particle i can have at any given time (left column) 
and update protocol that is followed according to the particle position after a simulation step (right columns). The arrays 
nR(i) and nL(i) store the indices of the disks that support disk i (see text). In case (d-4) we draw a small disk j since this 
situation is difficult to appreciate by eye for equal-sized disks where small differences in disk-to-wall overlaps are responsible 
for particle j moving to the right of i. 
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(3) either nR(i) or nL(i) is no longer in contact with i — > (i) set the lost contact to —2; (ii) if the remaining contact 
is nR(i) (the wall), then set nR(i) — —2, else check that nL(i) is still to the left of i and with a y-coordinate below 
particle i; if not, set also this contact to —2. 

(4) nR{i) and nL(i) are still in contact with i, but nL(i) is not to the left of i — > set nR(i) = nL(i) and set 
nL(i) = -2. 

(5) nR(i) and nL{i) are still in contact with i, nL(i) is still to the left of i, but the cm of i is not above the contact 
chord — ► set nR(i) and nX(i) to —2. 

Case (e,): This case corresponds to a particle i that was in the air (without defined supporting contacts) before the 
position update. We need to check if any contact has occurred and if it is a potential first stabilizing contact. Before 
doing this, we check that the particle has negative vertical velocity to avoid considering particles that are in their way 
up after a bounce. Two situations may arise: 

(1) Particle i contacts the container's base — ► set nR(i) and nL{i) to —1. 

(2) Particle i contacts another particle j — > if i has higher y-coordinate than j, then set either nR(i) or nL(i) 
(depending on which side of i is j) to j . 

Case (f): This case never arises because a contacting wall cannot be considered a first stabilizing contact. 

IV. RESULTS ON TAPPED GRANULAR BEDS 

In order to test the algorithm for the identification of supporting contacts and arches in a realistic non-sequential 
deposition we carry out MD simulations of the tapping of 512 disks in a rectangular box 13.91d wide. The values 
used for the parameters of the force model are: (i = 0.5, k n = 10 5 , 7 n = 300, k s = |fe n and 7 S = 200 with an 
integration time step S = 10 . The stiffness constants k are measured in units of mg/d, the damping constants 7 
in myj g/d and time in yj d/g. Here, m, d and g stand, respectively, for the mass of the disks, the diameter of the 
disks and the acceleration of gravity. The walls and base are represented by disks with infinite mass and diameter. 
The velocity- Verlet method was used to integrate the equations of motion along with a neighbor list to speed up the 
simulation. In order to keep stable the numerical integration, a basic requirement is to use an integration time step S 
much smaller than the contact time between particles, which is essentially controlled by 7„. We have chosen 5 to be 
50 times shorter than the typical duration of a contact. A more detailed discussion about the numerical algorithm 
can be found in Ref. [l4j . 

The value of 7„ is chosen deliberately high in order to have a small coefficient of normal restitution (e n = 0.058). 
This way we reduce computer time since fast energy dissipation leads the system to rest in less time. Of course, this 
also reduces the number of bounces and oscillations in the system. 

Disks are initially placed at random without overlaps in the simulation box, and the initial velocities are assigned 
from a Gaussian distribution of mean zero. The y-coordinates range from up to 100<i. Then, particles are allowed 
to deposit under the action of gravity. The deposit is said to be stable when the supporting contacts arrays nR(i) 
and nL(i) (see Sec. Ill), remain unchanged for 10 4 time-steps. Then, the stable configuration is saved. It is worth 
mentioning that this criterion to decide when the system is "at rest" is very reliable. All our simulation runs end 
up with unchanging supporting contacts arrays. Even though particles may perform small vibrations about their 
equilibrium positions and orientations, they remain supported by the same set of particles and this is what defines 
the mechanical stability in a macroscopic sense. 

The tapping process is simulated by moving the container's base and walls in the vertical direction half period of 
a sine function of amplitude A and frequency w = 2yJ g/d. The tapping amplitude is measured through the peak 
acceleration T — Aw 2 . After all particles have settled down and the stable configuration has been saved, we increase 
the amplitude of the perturbation by Ar = 0.0134y and the process is repeated. The range of the tapping amplitude 
goes from to 6.41y. When the maximum amplitude is reached the system is then tapped with decreasing amplitudes 
down to r = 0. Then, a new increasing T cycle is started and so on. 

The "annealing" protocol is performed several times starting from different initially deposited beds and then the 
results are averaged separating the very first increasing T phase from the rest of the tapping ramps. We introduce 
this averaging method due to the fact that all quantities obtained for the initial increasing ramp do not match the 
subsequent decreasing and increasing ramps for amplitudes below some threshold Tq [See inset of Fig. 3(b)]. The first 
increasing ramp of the tapping protocol is known as the irreversible branch and the subsequent coinciding ramps 
are known as the reversible branch. Despite the introduction of static friction (which could induce history dependent 
effects) and the fact that we apply a single tap to the system at each T, the results show the same trends as many 
experiments pH , ITtI Il8| where the granular assembly is subjected to some sort of tapped annealing. Unfortunately, 
the slow decaying time of the compaction process makes virtually impossible to simulate a real experimental situation. 

In order to enhance the quality of the averaged arch properties at given states along the tapping ramp we take 
intermediate configurations (r = 0.71 and 4.99) and tap the system at constant T for longer runs (1000 cycles). These 
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FIG. 2: Two snapshots of the configurations obtained for tapping amplitudes F — 0.71 (a) and 4.99 (b). Arches identified with 
the protocol of Sec. Ill are indicated by segments that join disks belonging to the same arch. Note that some disks may seem 
to be in unstable positions since all contacts are above its center. These disks are in fact stable thanks to the static friction 
forces. For each arch there are two disks that form the base of the arch, these are not indicated in the figure. 
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FIG. 3: Coordination number, support number and packing fraction as a function of the tapping amplitude F within the 
reversible branch of the annealing protocol, (a) Coordination number (z) (dashed line) and support number (z) Bupport (solid 
line), (b) Packing fraction cj>. The inset in panel (b) shows separately the curves corresponding to several tapping ramps 
including the irreversible branch (lower curve). 



values of V correspond to a high density and a low density state of the system. The results obtained from these longer 
runs are essentially the same as those obtained from the corresponding states along the annealing process; only the 
statistical dispersion is improved. 

We can see two snapshots of configurations obtained at low and high tapping amplitudes in Fig. |3 The arches 
identified in the sample can be appreciated by means of the segments that join mutually stabilizing disks (see Sec. 
III). A single disconnected segment indicates that the two joined disks support each other and therefore form an arch. 
A chain of connected segments indicates that all joined disks belong to the same arch. For each arch there are two 
disks that form the base of the arch, these are not indicated in the figure. The configuration with T = 0.71 [Fig.^a)] 
is quite ordered, showing a localized disordered region at middle height in the packing. Large crystal-like domains 
of triangular order are observed with clear defect boundaries in agreement with experiments 20]. Here, most arches 
consist of only two particles; only a few three-particle arches can be appreciated. Figure[2Ib) corresponds to T = 4.99. 
This configuration is clearly more disordered than the former; however, void spaces are rather uniform in size and 
homogeneously distributed throughout the packing. This contrasts with the larger degree of disorder displayed by 
pseudo dynamic simulations Q where a wider distribution of void sizes is found. Figure |2b) shows a larger number 
of arches than Fig. |2Ia), particularly larger arches consisting of up to six disks. 

In Fig.[3]we show the results of the mean coordination number {z) and packing fraction (f> obtained for the reversible 
branch of the tapping ramp (i.e., decreasing ramp). We also show the values for (z) support that includes only those 
contacts that serve to the stability of at least one of the two touching particles. Unlike hard particles Q, soft particles 
present more contacts than those just needed to make particles stable, therefore (2) support is not related to (z) in 
a trivial way. We have to bear in mind that is (z) support which is directly related to the arch size distribution and 
not (z). Our results show that (£) support and (z) present qualitatively the same behavior although (z) support presents 
lower values and varies within a narrower range. 

We can see from Fig. [3] that both (z) and <fi have rather high values when the tapping amplitude is small, corre- 
sponding to an ordered system. As T is increased those values undergo a slow decrease as the system gets disordered. 
Finally, <\> and (z) present a slight increase from their minimum values for high tapping amplitudes. This result qual- 
itatively agrees with the experimentally observed behavior in 3D 16] where a minimum in the density dependence 
on r has been observed. However, simulations of disks through pseudo dynamics show a much sharper transition 
between the ordered (low F) and the disordered (high T) regime. Moreover, these simulations |3| present a positive 
slope in (^) support for low tapping amplitudes within the ordered regime. 

Let us look into the relationship between the former quantities and the arches identified by our algorithm. Fig. 0] 
shows the number of arches normalized by the number of particles as a function of T along the reversible branch. 
We see that the large decrease in (z) and 4> (Fig. |3J) corresponds to an increase of the total number of arches in the 
bed. This number reaches a maximum value and then shows a slight decrease in correspondence with the increase 
in (z) and </> for high values of T. Clearly, there is a correlation between the buildup(brcakdown) of arches and the 
number of contacts and the free volume encountered in the sample. At very low tapping amplitudes the particles 
deposit in a very ordered manner without forming many arches; this means that, according to Eq. (Q, (z) t is 
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FIG. 4: Number of arches per disk as a function of F. Only the reversible (decreasing) branch of the tapping ramp is shown. 

rather high (ss 3.5, being 4 the maximum allowed value in 2D); moreover, there are few arch-trapped voids and 4> is 
high. At moderate T, arches are created with high probability; the coordination number then falls accordingly and 
the arch-trapped voids lower the packing fraction. For very high tapping amplitudes arches are again broken down 
partially; new particle-particle contacts are created and arch-trapped voids are filled up. 

In Fig. [S] the arch size distribution n(s) is shown for two values of the tapping amplitude. We can observe that 
for moderate tapping amplitudes the system has a larger number of large arches whereas for gentle tapping there is 
a larger amount of disks that do not belong to any arch (s = 1). Interestingly, in both cases the distribution can 
be fitted to a second order polynomial, in a semilogarithmic scale. This is in agreement with results from pseudo 
dynamics simulations 0. In 3D it has been found that n(s) oc s - 10 ± 003 which corresponds to a significantly 

higher prevalence of large arches with respect to 2D packings. 

In Fig. we show the distributions of the horizontal span of the arches and compare them with a theoretical model 
based on a restricted random walk. The horizontal span is the projection onto the horizontal axis of the segment 
that joins the centers of the right-end disk and the left-end disk of an arch. We use the results obtained for arches 
composed by two, three and four disks obtained for two different tapping amplitudes. At high packing fractions arch 
extensions appear discretized. Since the system presents an ordered layered structure, the particles of any arch of 
s disks form a string that connects layers (or stay in the same layer) so that the end-to-end horizontal projection 
of the arch can only be an integer number of on-layer jumps plus an integer number of between-layers jumps; both 
jumps have characteristic horizontal projections in the triangular lattice formed by the disks. We saw above that 
although <j) undergoes a decrease in its value for intermediate tapping amplitudes, it remains rather high. For this 
reason, although horizontal span distributions are more homogeneous, we still find some characteristic peaks for high 
r. These results do not agree with the prediction of the random walk model where horizontal span distributions are 
smooth. However, it must be taken into account that the restricted random walk model was proposed to represent 
arches at the outlet of a hopper where particles cannot order but simply fit in the gap between walls. 

V. ARCHES FROM STATIC CONFIGURATIONS 

In this section we intend to show to what extent an attempt to identify arches from the static structure of a deposited 
system can yield realistic results. We test two criteria that select two of the contacting particles of each particle as 
the supporting pair of the given grain: (a) random stabilizing pair, and (b) lowest stabilizing pair. We identify all 
pairs of contacting particles that may — according to their relative positions and the contact chord criterion — stabilize 
a given grain. For criterion (a) we then choose one of these pairs at random as the supporting pair. For criterion (b) 
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FIG. 5: Distribution of arch sizes for V — 4.99 (down triangles) and 0.71 (up triangles). The lines are fits to a second order 
polynomial. The data for s = 1 correspond to the number of disks that are not forming arches. 



we choose the pair that has the lower center of mass. In principle, we expect that case (b) should provide a more 
reliable identification of supporting contacts since those contacting particles in lower positions may correspond to the 
"real" supports, whereas those in higher positions may correspond to contacts made by further deposition of other 
particles on top of the particle under consideration. We will see that this is not the case. 

We have obtained the arch size distribution n(s) by applying the static criteria above to packings corresponding 
to r = 0.71 and 4.99. Fig. [7\ shows clearly that, for disordered packings (r = 4.99), arches are identified much 
better by the random criterion. The lowest supporting pair criterion detects less arches than the realistic dynamic 
criterion. As we can see from the inset of Fig. for very ordered packings (r = 0.71), both criteria fail. Criterion 
(a) overestimates whereas criterion (b) underestimates the number of arches beyond s — 2. A detailed analysis of 
the supporting contacts detected by the static criteria shows that criterion (a) tends to find many false mutually 
stabilizing pair of particles in ordered packings, which leads to the identification of false arches. On the other hand, 
criterion (b) tends to detect too few of the real mutually stabilizing pairs of particles and for that reason the number 
of arches detected are fewer than with the dynamic criterion. 



VI. CONCLUSIONS AND FINAL REMARKS 



We have presented a protocol to dynamically identify arches during the deposition of granular particles carried 
out through realistic granular dynamics. We find that simple static criteria can correctly identify arches to some 
degree for not very ordered packings. In particular, the criterion we call random stabilizing pair seems to be the 
most successful and could be used to identify arches in experimentally generated 2D granular beds |l9ll20| . However, 
in ordered packings is very difficult to identify arches from static configurations. The reason for this is that in an 
ordered packing the coordination number is rather high and then the number of plausible stabilizing pairs for a given 
particle is higher than that in disordered packings. Since choosing the correct pair is the essence of the criterion to 
successfully identify stabilizing contacts, the more pairs we have to chose from, the more likely we make a mistake. 

The packings generated by varying the tapping amplitude "quasistatically" present a low-density disordered regime 
(at high tapping amplitudes) and a high-density, ordered regime (at low tapping amplitudes). The support number 
decreases with increasing tapping amplitude except for of the very high tapping amplitudes. These observations 
are in contrast with pseudo dynamic simulations Q where (z) SU p por t increases with T (except at the order-disorder 
transition) . We believe that this is due to the fact that the pseudo dynamics is a good representation for fully inelastic 
disks that roll without slip and that "fall down" at constant velocity — a situation that seems to be met by particles 
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FIG. 6: Horizontal span distribution for arches with 2 (a), 3 (b) and 4 (c) disks. The lines correspond to V = 0.71 (dotted) 
and r = 4.99 (dashed). The solid lines correspond to the restricted random walk model Q 
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FIG. 7: Distribution of arch sizes for T = 4.99 comparing the "real" distribution obtained from simulations (circles) with the 
distributions obtained from the static structure using criteria (a) (up triangles) and (b) (down triangles). The inset shows the 
same results for T — 0.71. 



carried on a conveyor belt at low velocities [19( . The simulations in the present work are rather far from this regime 
since particles can spring away from an impact and so explore a wider range of position in the box before finding a 
locally stable configuration. 

The form of the arch size distributions is in agreement with pseudo dynamic simulations. This suggests that the 
form of n(s) is rather insensitive to the deposition algorithm. The calculated number of arches grows sharply as the 
density decreases when the tapping amplitude is increased. However, this tendency is not monotonic and above a 
transition zone the number of arches decreases whereas the density increases. This feature highlights the intuitive 
connection between arches and density fluctuations. 

We have to point out here that the criterion we used to decide if two contacts may be the stabilizing contacts for a 
given particle i (i.e. that the contact chord is below the cm of i), rules out the possibility that a particle be supported 
from "above". This situation does indeed arises occasionally due to static friction. Two contacting particles with 
y-coordinates slightly above particle i may sustain the particle due to high static friction forces. We have found that 
this situation indeed occurs in our simulations. Indeed, some instances can be appreciated in Fig. |^b). A detailed 
analysis of these types of arches that are also found in experiments (see for example Ref. Q) will be presented 
elsewhere. 
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